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PREFACE 


This  report  is  one  of  a series  documenting  NBS  research  and  analysis  efforts  to 
support  the  Department  of  Energy/National  Bureau  of  Standards'  Measurements 
Program.  The  work  reported  in  this  document  was  performed  cooperatively  and 
supported  by  DoE/NBS  Task  Order  A008  under  Interagency  Agreement  No.  DE-AIOl- 
76PR06010. 
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ABSTRACT 


A discussion  is  presented  of  theory  of  heat  transfer  in  low-density,  glass- 
fiber  insulation  via  conduction,  convection,  and  radiation.  It  is  concluded 
that  the  primary  modes  of  heat  transfer  in  this  material  are  air  conduction  and 
radiation.  An  analysis  of  NBS  data  of  measured  apparent  thermal  conductivity 
for  different  thicknesses  results  in  a parameter  estimate  of  the  optical 
extinction  coefficient.  This  parameter  determines  the  amount  of  change  in 
apparent  thermal  conductivity  as  a function  of  sample  thickness.  This  phenomena 
is  referred  to  as  the  "thickness  effect." 

Keywords:  conduction;  convection;  guarded-hot-plate;  low-density  insulation; 

thermal  conductivity;  thermal  radiation;  thickness  effect. 
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Nomenclature 


A = area  (m^) ; 

e = blackbody  emissive  power,  oT‘^(W»m~2); 

fy  = volume  fraction,  defined  by  equation  (3); 

g = gravitational  constant,  9.8  (m»s“2); 

= non-dimensional  temperature  integral  defined  in  equation  (25); 
K = permeability  (m^); 

k = thermal  conductivity  (W*m“^ ; 

L = thickness  of  insulation  (m); 

n = iteration  counter; 

Nu  = Nusselt  number; 

P = porosity; 

6 = heat  flow  (W); 

q = heat  flux  (W*m~2)j 

= ratio  of  conductivities,  k^/kg; 

Ra*  = modified  Rayleigh  number,  defined  by  equation  (5); 
r2  = goodness  of  fit  parameter,  defined  by  equation  (9); 

T = temperature  (K); 

AT  = temperature  difference,  T]^  - T2(K)  ; 

X = non-dimensional  spatial  coordinate,  X'/L; 

x’  = dimensional  spatial  coordinate  (m) ; 

= first  unknown  parameter  (kg»m  ); 

f 

Xj^  = linearized  (old)  value  of  the  first  parameter  (kg*m  ); 
y^  = measured  value; 

Greek  Symbols 

= stagnant  thermal  diffusivity  (m  *s  ); 

vi 


3 = mass  attenuation  coefficient  (m^'kg  or  coefficient  of  cubical 

expansion  (K~l); 

Y = second  unknown  parameter  (W*m"l*K~l); 

5 = deviation; 

e = emissivity; 

V = kinematic  viscosity  (m-*s“M; 

p = density  (kg*m~3); 


a = standard  deviation; 

Tq  = optical  depth,  defined  by  equation  (12). 

Subscripts 


a 

denotes 

air  or  "in  the  air" 

c 

denotes 

convection  or  critical 

e 

denotes 

effective 

g 

denotes 

glass 

i 

denotes 

insulation 

j 

denotes 

j * value 

k 

denotes 

conduction  or  k’^^^  iteration 

1 . s . 

denotes 

least  squares 

m 

denotes 

metering  section 

r 

denotes 

radiation 

t 

denotes 

total  (i.e.  the  entire  48"  x 48 

w 

denotes 

at  the  wall 

0 

denotes 

stagnant 

1 

denotes 

the  hot  surface 

2 

denotes 

the  cold  surface 

sample) 


vii 


Table  of  Contents 


Page 

PREFACE  ill 

ABSTRACT  iv 

ACKNOWLEDGMENT  v 

NOMENCLATURE  vi 

1.  THICKNESS  EFFECT  IN  LOW-DEJISITY  INSULATION  1 

1.1  INTRODUCTION  1 

2.  HEAT  TRANSFER  OVERVIEW  4 

2.1  THEORETICAL  FORMULATION  4 

2.2  CONDUCTION  4 

2.3  CONVECTION  5 

2.4  RADIATION  7 

3.  ANALYSIS  9 

3.1  PARAMETER  ESTIMATION  9 

3.2  SEOUENTIAL  FORMULATION  10 

3.3  APPLICATIONS  12 

4.  CONCLUSIONS  AND  RECOMMENDATIONS  19 

5.  REFERENCES  21 

APPENDIX  A A-1 


viil 


1.  THICKNESS  EFFECT  IN  LOW-DENSITY  INSULATION 


1.1  INTRODUCTION 

Determining  the  heat  transfer  through  porous  insulating  materials  has  long  been 
a subject  of  concern.  In  particular,  judicious  engineering  design  requires  the 
minimization  of  the  total  cost  associated  with  the  purchase  of  insulation. 

This  total  cost  includes  the  base  cost  of  the  insulation  and  the  total  cost  of 
the  energy  that  flows  through  the  insulation.  The  base  cost  consists  of  the 
material  cost  and  the  installation  cost.  Furthermore,  the  cost  of  the  energy 
that  is  transmitted  through  the  insulation  must  be  amortized  over  the  lifetime 
of  the  insulation,  and  it  is  convenient  to  perform  this  minimization  in  terms 
of  net  present  worth. 

Indeed,  to  effect  the  minimization,  it  is  necessary  to  characterize  the 
performance  of  the  insulation  with  parameters  that  can  be  associated  with  cost. 
Thus,  the  thermal  performance  of  insulating  materials  is  presently  characterized 
versus  the  density  of  the  sample  and  the  mean  diameter  of  the  fibers  comprising 
the  porous  matrix  [1],  Moreover,  such  has  been  the  practice  for  some  time,  as 
many  have  analytically  addressed  [2-5]  the  problem  of  modeling  the  thermal  per- 
formance of  insulating  materials,  and  the  experimental  investigations  are 
legion  [4-8]  also.  Recently,  the  research  emphasis  has  been  dedicated  towards 
determining  the  "thickness  effect"  [8-11],  that  is,  the  change  in  the  effective 
thermal  conductivity  with  specimen  thickness.  In  particular,  a very  detailed 
discussion  of  the  thickness  effect  may  be  found  in  Reference  11.  In  connection 
with  this,  it  is  noteworthy  that  the  effective  thermal  conductivity  can  be  well- 
characterized  [12]  versus  percent  light  transmission. 
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To  date,  the  problem  has  been  addressed  only  deterministically  in  that  physical 
mechanisms  have  been  assumed,  ranges  of  transfer  coefficients  have  been 
selected,  and  solutions  have  been  obtained  under  a wide  variety  of  simplifying 
assumptions.  This  technique  yields  excellent  insight  when  determining  the 
relative  effects  of  various  physical  mechanisms,  but  it  is  sometimes  fruitless 
to  attempt  to  accurately  predict  the  actual  heat  transfer  occurring  in  a speci- 
fic physical  system  because  the  predicted  property  (in  this  case  the  effective 
thermal  conductivity)  may  vary  much  more  under  the  range  of  input  parameters 
than  does  the  observed  experimental  data.  Actually,  this  can  be  attributed 
both  to  the  variation  due  to  the  approximate  nature  of  the  solutions,  and  due 
to  the  lack  of  knowledge  of  the  specific  input  parameters.  In  particular,  the 
solid  conduction  contribution  will  change  depending  upon  the  assumptions  made 
as  to  fiber  orientation,  and  the  radiative  parameters  (phase  function  and 
extinction  and  single  scattering  coefficients)  are  not  well  known  for  general 
situations.  The  adhesive  binder  is  a complicating  factor  of  great  physical 
importance  that  is  usually  neglected  in  theoretical  analyses,  and  there  is 
likely  to  be  large  spatial  variations  in  binder  content  in  any  one  batch  of 
insulation.  There  is  also  likely  to  be  large  global  variations  in  binder  con- 
tent from  batch  to  batch.  Thus,  the  law  of  diminishing  returns  is  operating 
in  that  further  deterministic  analyses  will  be  very  tedious  and  will  yield 
less  practical  information  than  those  of  the  past. 

In  this  light,  a stochastic  approach  will  be  undertaken  herein  that  will 
circumvent  many  of  these  difficulties  and  yield  some  new  insight  into  the 
experimental  determination  of  the  effective  thermal  conductivity.  In  short, 
relationships  for  heat  transfer  modes  will  be  assumed  in  terms  of  unknown 
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constants.  These  unknown  constants  will,  in  effect,  be  the  uncertain  parameters 
mentioned  above.  Since  these  parameters  will  be  determined  from  the  data,  they 
will  automatically  be  the  "best-fit”  values  consistent  with  the  assumed  model 
and  the  scatter  in  the  data.  The  true  test  of  this  method  will  entail  compari- 
son of  the  overall  scatter  in  the  data  against  the  scatter  in  the  predicted 
values.  Also,  enough  data  points  must  be  used  to  ensure  a statistically  signi- 
ficant result.  Both  of  these  facets  of  the  present  approach  will  be  quantified. 

Because  the  parameters  will  be  derived  from  experimental  data  obtained  at  the 
NBS  in  Gaithersburg,  MD,  a few  comments  concerning  the  quality  and  the  material 
variability  of  the  insulation  samples  are  in  order.  Figure  1,  which  was  abstracted 
from  Reference  13,  shows  that  both  density  and  fiber  diameter  are  important 
parameters  in  determining  the  effective  thermal  conductivity  of  fibrous  insulations. 
The  data  used  herein  are  NBS-measured  values  taken  from  the  calibrated  transfer 
samples  prepared  for  the  NBS  by  the  Johns  Manville  Corporation.  These  samples 
vary  widely  in  thickness  and  density,  yet  the  fiber  diameters  were  controlled 
very  carefully  to  have  a mean  value  of  3.8  x 10“^  m for  all  of  the  samples. 

Further  detail  may  be  found  in  Reference  1.  Thus,  we  shall  proceed  realizing 
that  the  fiber  diameter  is  a parameter  which  is  held  fixed  in  the  present 
analysis.  The  most  that  can  be  done,  then,  is  to  predict  the  thickness  effect 
at  a constant  fiber  diameter.  For  completeness.  Table  1 includes  values  of 
the  datum  points  used. 

The  present  analysis  will  proceed  in  the  following  vein.  Physical  mechanisms 
for  the  heat  transfer  will  be  examined  and  proposed,  and  detailed  discussions 
of  past  models  are  given.  A brief  description  of  the  sequential  least  squares 
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algorithm  follows.  The  data  are  then  filtered  through  a variety  of  models, 
and  the  present  work  ends  with  conclusions  and  recommendations. 

2.  HEAT  TRANSFER  OVERVIEW 

2.1  THEORETICAL  FORMULATION 

The  problem  under  consideration  is  to  predict  the  heat  flow  through  low  density 
insulation  given  the  boundary  conditions  (the  boundary  temperatures).  Of 
course,  the  basic  material  properties  must  also  be  known  as  they  are  required 
input  parameters.  The  heat  flow  will  be  assumed  to  be  one-dimensional  as  the 
NBS  GHP  (Guarded  Ifot  Plate)  satisfies  very  stringent  criteria  in  this  regard 
(about  0.1  percent  accuracy).  In  this  case,  the  heat  flow  may  be  conceptualized 
as  follows; 

Q = + Qc  + 0^  0]^c  ^r  ^c 

where  separate  contributions  have  been  included  for  conduction  (k) , convection 
(c),  and  radiation  (r).  Also,  interactions  between  all  three  modes  of  heat 
transfer  have  been  allowed  as  separate  terms  and  it  is  presumed  that  these  terms 
are  of  smaller  magnitude  than  the  ones  involving  the  basic  modes  of  heat 
transfer.  The  validity  of  this  assumption  will  be  demonstrated  by  the  match  between 
thing  and  experiment.  The  terms  will  be  addressed  one-by-one. 

2.2  CONDUCTION 

The  conduction  terra,  Q[^,  involves  the  gaseous  and  solid  conduction  through  the 
porous  medium.  There  is  also  an  interaction  term  associated  with  these  two 
modes  of  conduction.  This  can  be  conceptualized  as  follows:  since  the  contact 

areas  between  fibers  are  small,  the  constriction  resistance  is  large  and 
little  heat  flows  through  the  solid  matrix  directly.  However,  if  the  fibers  do 
contact  or  if  they  come  close  together,  heat  will  flow  from  one  fiber  through  the 
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intervening  gas  to  another  fiber.  A careful  analysis  of  this  process  for  packed 
spheres  [14]  shows  that  this,  in  fact,  is  a major  transfer  mechanism  for  heat 
flowing  through  the  solid.  Indeed,  a number  of  analytical  investigators  [2-5, 

9,  10,  14-16]  have  pursued  this  subject.  Experimental  treatments  of  this  topic 
[2-8]  confirm  most  of  the  theories  for  low  density  porous  insulation  wherein 
the  solid  matrix  is  so  widely  dispersed  (porosity  = 0.995)  that  only  a small  por- 
tion of  heat  flow  can  be  attributed  to  the  fibers.  This  is  readily  evident  by 
viewing  Figure  2,  which  has  been  abstracted  from  the  review  article  by  Bankvall 
[5].  The  model  used  herein  is  one  proposed  by  Bhattacharyya  [10]: 


1 -_ 


1 - P 


1 4-  2 f^ 
1 + Rk 


+ (l-\)  [1  - 


1 - R, 


(2) 


1 + fv(  1 + 5R^) 

3(1  + Rk) 


where,  for  air  at  287K  and  glass  fibers,  the  properties  are  taken  as 
follows:  kg  = 1.04  Wm~^‘K~^,  k^  = 0.026  W*m^  • K“^  , and  Pg  = 2580  kg*m"^, 
p_  = 1.04  kg*m  . The  volume  fraction  is  computed  as  follows: 
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2.3  CONVECTION 


f = 


(3) 


The  convective  heat  transfer  in  horizontal  porous  media  has  been  the  subject  of 
extensive  investigation  for  several  decades  [17-30].  The  critical  modified 

•k 

Rayleigh  number,  Ra^ , below  which  no  convection  can  occur  in  a horizontal 
porous  medium  heated  from  below  has  been  determined  to  be  as  follows  [17,  18]: 
Ra*  = 4r2  2 39.5  (4) 

where  the  modified  Rayleigh  number  qualitatively  represents  the  ratio  of  the 
buoyancy  force  to  the  drag  force  of  the  porous  matrix,  and  is  given  as: 
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(5) 


= g g at  K L 
vao 

wherein  the  effective  thermal  conductivity  of  the  stagnant  porous  medium  is 
used  to  calculate  the  thermal  diffusivity,  o.q,  A number  of  experimental 
investigations  [20,  22-24]  have  verified  the  relationship  given  by  equation  (4). 
Elder  [24],  gives  the  following  formula  for  the  non-dimensional  steady  convec- 
tive heat  transfer  in  porous  media  heated  from  below: 

• • 

Q + Q * * 

Nu  = S Ji  = 1 + (Ra  - Ra(,)/40  (6) 

where  Nu  = 1 if  Ra  ^ ^c* 

It  is  of  interest  to  estimate  the  magnitude  of  the  free  convective  heat  transfer 
in  porous  media  using  equations  (4)  and  (5).  First  the  modified  Rayleigh  num- 
ber is  calculated.  The  following  values  of  the  parameters  are  used: 
g = 9.8  m-s“^,  8 = T“^  = 3.37  x 10“^  K“^,  v = 1.56  x lO"^  m^*s"^ 
k*  = 0.045  W*m“^*K“^,  p = 1.16  kg»m“^,  Cp  = 1.011  J*kg“^»K“^,  = 

3.84  X 10“5  L = 0.1562  m,  K = 5 x 10~®  m^.** 

These  values  result  in  a modified  Rayleigh  number  of: 

Ra*  (AT)  = 0.419  AT(K)"1  (7) 

For  the  AT  = 27. 8K  used  in  the  GHP,  the  modified  Rayleigh  number  is: 

Ra*  (27.8)  = 11.6  (8) 

which  is  below  the  critical  value,  so  all  of  the  terms  associated  with 
convection  in  equation  (1)  are  nil.  Indeed,  setting  equation  (7)  equal  to 

* Taken  from  measurements  made  using  the  NBS  Guarded  Ifot  Plated  [1]. 

**  Taken  from  Fournier  and  Klarsfeld  [7],  Fig.  5,  F = 3.6  (typical  of  the 

3.8  pm  diameter  fibers  in  the  measurement  system  of  [1]). 
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equation  (4)  results  in  a AT  = 95. 9K  for  the  onset  of  convective  flow.  Thus, 
free  convection  is  seen  to  be  a significant  contributor  to  the  heat  transfer 
for  porous  insulations  only  at  very  large  thicknesses  and  very  large  temperature 
differences.  Indeed,  for  all  practical  cases  of  thermal  insulation  at  room 
temperature,  free  convection  heat  transfer  does  not  occur  in  the  horizontal 
geometry . 

2.4  RADIATION 

As  seen  by  viewing  Figure  2,  radiation  heat  transfer  contributes  significantly 
to  the  total  heat  transfer  in  low  density  fibrous  insulation.  Recent  theoreti- 
cal calculations  neglecting  the  phenolic  binder  [31]  show  that  the  scattering 
coefficient  is  roughly  a factor  of  2 to  4 times  larger  than  the  absorption 
coefficient,  but  this  conclusion  has  not  been  conclusively  supported  with 
experimental  data.  Indeed,  there  is  some  indication  that  the  phenolic  binder 
is  highly  absorptive.  Moreover,  it  is  extremely  difficult  to  measure  absorption 
and  scattering  coefficients  separately  [32-33],  especially  in  the  midst  of  a 
matrix  of  solid  fibers.  Despite  these  limitations,  a brief  review  of  past 
work  in  this  area  will  be  presented  for  illustrative  purposes. 

The  pioneering  work  was  done  by  Hamaker  [34]  and  Larkin  and  Churchill  [35]  who 
applied  the  two  flux  model  for  gray  absorbing  and  scattering  media.  The  approx- 
imation of  radiative  equilibrium  used  does  not  allow  free  interaction  between 
the  radiation  and  conduction  which  limits  the  general  applicability  of  their 
results.  Viskanta  subsequently  pursued  the  problem  including  isotropic  scatter- 
ing from  a more  rigorous  standpoint  for  the  condition  of  radiative  equilibrium 
[36],  and  later  coupled  with  conduction  [37],  It  is  mentioned  in  passing  that 
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Usiskin  and  Sparrow  [38],  and  Heaslet  and  Warming  [39]  presented  a solution  for 
a purely  absorbing  medium.  Furthermore,  the  diffusion  technique  has  also  been 
applied  to  the  problem  [46]  of  radiation  heat  transfer  in  a purely  absorbing 
medium. 

Some  further  analyses  into  the  problem  at  hand  have  attempted  to  remove  the 
assumption  of  isotropic  scattering  [41-43],  enhance  the  applicability  of  the 
two-flux  assumption  [44,  45]  or  include  transients  [46].  Further  analytical 
work  may  be  found  in  Reference  47,  wherein  a large  parameter  variation  was 
performed.  All  of  the  investigations  mentioned  so  far  assume  gray  material 
properties . 

Recently,  Tong  and  Tien  [31]  have  applied  a two-flux  model  to  the  specific 
problem  involving  glass  fibers  wherein  the  spectral  solution  of  Maxwell's  equa- 
tions have  been  used  to  determine  the  absorption  and  scattering  coefficients. 
While  this  represents  a large  step  in  the  right  direction,  the  questions  of 
fiber  orientation  and  the  effect  of  the  adhesive  binder  are  still  to  be  included 
in  a more  rigorous  analysis.  Furthermore,  it  is  questionable  whether  such  an 
approach  involving  an  approximate  solution  technique  can  ever  be  predictive  to 
within,  say,  two  percent  as  is  necessary  to  predict  the  thickness  effect. 

It  is  therefore  concluded  that  the  research  community  has  yet  to  accurately 
solve  the  generalized  equation  of  transfer  for  the  specific  problem  of  radiant 
heat  transfer  in  fibrous  insulations  with  binder  present.  Furthermore,  the 
property  measurements  that  have  been  made  [e.g.,  extinction  coefficient]  must 
be  augmented  by  further  measurements  including  the  single  scattering  albedo 
[32-33]  to  be  generally  applicable  to  the  present  problem.  Indeed,  this  is 
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not  an  easy  problem  either  as  the  question  of  the  separation  of  effects  when 
measuring  quantities  that  can  be  attributed  to  diverse  physical  mechanisms  is 
an  elusive  one  [48]. 

Rather,  the  approach  adopted  here  will  consist  of  taking  the  measured  data  for 
heat  transfer  rates,  assuming  physical  models  and  filtering  the  measured  data 
through  the  assumed  model  to  attain  "best-fit"  values  for  the  thermal  parameters. 
This  approach  has  the  advantages  of  yielding  the  best  thermal  parameters  that 
pertain  to  the  data  at  hand  consistent  with  the  noise  level  of  the  data,  the 
assumed  model  and,  of  course,  the  minimization  of  the  L-2  norm.  However,  the 
approach  cannot  truly  be  called  deterministic  as  the  material  properties  are 
derived  from  performance  data  rather  than  vice  versa,  but  this  is  a philosophical 
and  not  a technical  issue  and,  consequently,  will  not  be  further  addressed. 
Instead,  the  pitfalls  inherent  in  such  an  attempt  will  be  discussed,  and  it 
will  become  evident  that  the  range  of  applicability  of  the  generated  parameters 
will  be  pronouncedly  affected  by  the  specific  models  assumed  and  by  the  specific 
data  filtered. 

3.  ANALYSIS 

3.1  PARAMETER  ESTIMATION 

When  parameter  estimation  is  employed  to  determine  thermal  parameters,  two 
general  rules  must  be  separately  observed  in  order  to  obtain  meaningful  results. 
Firstly,  the  model  must  be  representative  of  the  physical  processes  occurring 
in  the  experimental  system.  For  example,  problems  are  sometimes  encountered 
when  free  convection  (T^/^  dependence)  and  radiation  (T^  dependence)  are  erron- 
eously lumped  into  the  same  heat  transfer  coefficients  in  a model.  This  type 
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of  model  will  not  allow  modes  to  be  identified.  Also,  all  terms  of  signifi- 
cance must  be  included  in  the  model.  If  convection  is  indeed  an  important 
mode  of  heat  transfer,  its  effects  will  be  contained  in  the  data,  and  when  the 
data  are  analyzed,  these  effects  will  be  attributed  to  some  parameters,  and 
incorrectly  so  if  convection  is  not  represented  in  the  model. 

Secondly,  the  combination  of  the  noise  level  of  the  data  and  the  number  of  data 
samples  taken  must  be  such  as  to  render  the  results  statistically  Tiomogeneous . 

The  general  rule  of  thumb  is  that  the  noise  level  of  the  data  should  be  much  less 
than  the  signal  level  of  the  physical  mechanism  that  is  to  be  identified.  This 
stringent  criterion  may  be  relaxed  as  larger  numbers  of  datum  points  become 
available.  Other  problems  that  are  specific  to  each  individual  problem  under 
investigation  [47]  will  not  be  addressed  here. 

3.2  SEQUENTIAL  FORMULATION 

The  formulation  of  the  sequential  linear  least  squares  algorithm  is  well-developed 
elsewhere  [48]  and  will  not  be  repeated  here.  Instead,  the  use  of  the  technique 
will  be  presented  by  example.  Suppose  that  j parameters  xj  are  to  be  identified 
from  the  measured  values  A,  which  multiply  the  parameters  and  y,  which  appear 
separately.  Mathematically, 

Akj  Xj  = yk  (9) 

where  k denotes  the  k’th  measurement.  The  system  becomes  overdeterrained  when  k 
becomes  greater  than  j,  and  no  general  solution  is  then  possible.  Rather,  the 
system  must  be  solved  for  the  parameters  that  satisfy  some  global  criterion. 

It  is  convenient  to  choose  this  criterion  to  be  the  minimization  of  the  L-2 
norm  of  the  errors  in  the  parameters. 


10 


When  identifying  j parameters,  the  linear  system  of  equations  becomes  complete 
when  j experimental  observations  have  been  added  to  the  system.  The  equations 
can  be  solved  at  this  time  by  Gaussian  elimination,  yielding  a set  of  parameters 
that  are  representative  of  the  measurements  taken.  When  these  measurements 
contain  significant  noise,  the  parameters  are  likely  to  be  substantially  in 
error  at  this  point.  In  order  to  yield  physically  realistic  parameters, 
additional  data  points  are  processed. 

However,  as  one  more  observation  is  added  to  the  set,  the  system  becomes 
overdetermined.  Then,  least  squares  procedures  are  used.  Least  squares  data 
processing  will  yield  a new  set  of  parameters  each  time  an  additional  observation 
is  processed.  This  is  possible  to  do  with  ordinary  least  squares  algorithms  by 
reformulating  and  reducing  the  A matrix  each  time.  The  sequential  least  squares 
algorithm  performs  this  function  at  a fraction  of  the  computer  time  because  the 
matrix  does  not  get  reformulated  and  reduced  for  each  data  point.  Thus, 
sequential  least  squares  yields  values  for  the  parameters  after  each  data 
observation  is  added. 

This  is  advantageous  in  two  ways.  First,  by  observing  the  parameter  trajectories, 
it  can  usually  be  determined  whether  enough  data  points  have  been  processed  for 
the  noise  to  be  filtered  out.  Secondly,  it  can  sometimes  be  determined  whether 
there  are  imperfections  in  the  model.  To  explain  this,  consider  the  situation 
where  data  have  been  taken  on  a 0.05  m specimen.  Suppose  that  the  initial  noise 
in  the  parameters  has  died  out,  and  the  parameter  trajectories  have  been  observed 
to  become  flat.  Suppose,  further,  that  data  from  a 0.1  m sample  are  next  added. 

If  the  parameter  trajectories  are  observed  to  suddenly  start  to  change  and 
continue  to  change  steadily,  it  can  be  deduced  that  there  is  error  in  the  model 
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as  it  does  not  fit  all  the  different  physical  situations  in  a similar  fashion. 
This  is  termed  "model  noise,"  Thus,  the  model  should  be  reformulated  and  the 
data  reprocessed  until  a uniformily  valid  set  of  parameters  is  obtained. 

In  summary,  a sequential  observation  of  the  parameter  trajectories  can  be  useful 
in  detecting  model  noise  and  in  determining  when  a sufficient  number  of  data 
observations  have  been  processed  to  overcome  the  experimental  noise  level.  The 
sequential  algorithm  is  given  for  reference  in  the  Appendix, 

3,3  APPLICATIONS 

The  problem  at  hand  is  formulated  in  terms  of  the  effective  thermal  conductivity 
as  follows: 

^e  ~ ^r  (10) 

where  is  the  term  in  brackets  in  equation  (2),  It  is  for  pure  conduction 
thru  the  solid  and  air. 

The  radiation  conductivity  is  taken  to  be  as  follows: 

2 2 

a(T,  + T-)(T,  + T-)L 

kj.  = i ± ± (11) 

L+L  -1+3t/4 
£l  S2 

where  all  of  the  quantities  in  equation  (11)  are  assumed  to  be  measured  except 
for  Tq,  the  effective  optical  depth.  This  parameter  is  termed  effective  in  that 
it  incorporates  the  effects  of  absorption,  anisotropic  scattering,  spectral 
properties,  fiber  orientation  and  the  binder  material.  It  is  standard  practice 
to  determine  the  properties  of  fibrous  insulation  (i.e,  the  effective  thermal 
conductivity  and  the  permeability)  as  a function  of  the  material  density,  p. 
Since  the  scattering  of  radiation  is  proportional  to  the  number  of  scattering 
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particles  within  the  medium,  it  is  reasonable  to  assume  that: 


Tq  = 3 pL  (12) 

Further,  it  is  known  [1]  that: 

L + - 1 «±^o  (13) 

£2  ^ 

Equations  (11)  and  (12)  may  be  combined  together  with  the  condition  stipulated 
as  equation  (13)  and  the  binomial  theorem  to  yield: 


4-(X?  + T?)(Ti  + T^)  . 4(1/£^  + 1/e2  - 1). 

- >^1^°  3P  ^ J 


where  the  term  has  been  linearized  by  the  inclusion  of  the  old  value  of  the 


parameter  Also: 


(14) 


= 3"^  (15) 

It  is  apparent  that  it  is  necessary  to  iterate  to  determine  a "thickness  effect" 

f 

as  it  only  shows  up  in  equation  (14)  as  a product  involving  the  old  value,  x^, 

f 

of  the  radiation  parameter.  Thus,  is  set  to  zero  (the  first  old  value), 

f 

Xj^,  is  determined;  the  old  value  x^^  is  next  set  to  this  value  of  Xj^ ; the 
filtering  is  performed  again  yielding  a new  value  of  xi  and  the  procedure 
continues  until  convergence  is  achieved.  Finally,  it  is  noted  that  equation 
(11)  is  similar  in  form  to  the  derivations  of  Verschoor  and  Greebler  [2], 
Bankvall  [5],  van  der  Held  [50],  Hager  and  Steere  [51],  and  Strong  et.  al . 


[52], 


The  geometry  for  the  first  problem  is  shown  in  Figure  3.  At  the  low  densities 
at  which  the  present  data  exist,  solid  conduction  plays  very  little  role  in 
the  heat  transfer,  and,  therefore,  the  conductive  contribution  was  modeled  as 
in  equation  (2),  and  taken  as  a known  quantity  (i.e.,  included  in  the  right- 
hand  side)  . 
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The  linearized  multiplier  for  the  radiative  parameter  is  explicitly: 


4o(Ti  + i)(Ti  + i2)  r,  - I’l 

Akl  3^  11  1=1  J 

Finally,  the  known  quantity  including  the  effective  thermal  conductivity, 
is  explicitly: 

yk  = ke  ~ ^ L/(AT  A) 

where  the  second  term  (=  k]^)  is  calculated  using  equation  (2). 


(16) 
ke  > 


(17) 


Figure  4 shows  the  variation  of  X]^  , with  the  number  of  datum  points,  n.  Recall 
that  a sequential  filter  generates  the  "best-fit"  value  for  the  parameter  being 
estimated  for  each  sequence  of  datum  points  up  to  and  including  the  present  one. 
The  data  are  organized  from  left  to  right  as  the  0.0254  m (1")  set,  the  0.0381  m 
(1.5")  set,  the  0.0762  m (3")  set  and  the  0.1524  m (6")  set.  Local  variations 
are  seen  to  be  extreme  below  n = 10  and  are  observed  to  damp  out  for  n>20  due 
to  the  averaging  of  a larger  number  of  points.  Thus,  it  may  be  concluded  that 
there  are  more  than  enough  data  points  to  render  local  fluctuations  (measurement 
noise)  insignificant,  However,  the  global  variation  (model  noise)  of  the 
parameter  manifested  as  a decrease  in  x^  with  an  increase  in  thickness  (negative 
slope)  is  disturbing.  This  suggests  that  the  proposed  model  for  the  effective 
thermal  conductivity  is  of  the  wrong  form  else  the  parameter,  X]^ , would  have 
leveled  off.  Obviously,  relative  to  the  value  of  x^,  at  a thickness  of  0.0254 
m (1"),  the  value  of  x^,  for  a thickness  of  0.1524  m (6")  should  be  lower. 

Thus,  although  the  family  of  predicted  curves  should  nestle  close  to  the  average 
of  the  experimental  data,  the  thickness  effect  will  not  be  well  represented. 


14 


This  is  indeed  observed  in  Figure  5 which  presents  the  predicted  curves  and  the 
experimental  data  points  for  the  four  thicknesses  of  interest.  While  the  aver- 
age values  have  been  predicted  well,  the  thickness  effect  has  been  overpredicted. 
The  0.0254  m (1")  curve  lies  on  or  below  all  of  the  0.0254  m (1")  data  but  one 
point.  The  0.0381  m (1.5")  curve  is  statistically  unimportant  as  there  are 
only  3 data  points  appertaining  thereto.  The  0.0762  m (3")  curve  appears  to 
be  a good  fit  for  its  data  and  the  0.1524  m (6")  curve  lies  above  all  of  its 
associated  data  points  but  one.  For  the  data  shown,  the  thickness  effect  is 
eyeballed  to  be  about  3.5  percent  at  p = 9kg»m“3)  while  the  predicted 
curves  indicate  a thickness  effect  of  about  6 percent. 

Several  comments  are  appropriate  here.  First  of  all,  the  value  of  3 = 39.5 
m2*kg“l  is  the  best  value  of  3 consistent  with  the  present  model  for  the  given 
experimental  data.  Secondly,  even  though  the  scatter  in  the  data  is  substantial, 
it  is  apparent  from  Figure  6 that  the  thickness  effect  is  overpredicted. 

Finally,  the  additional  correction  required  to  bring  the  predictions  in  line 
with  the  data  is  small  (about  2.5  percent),  and,  therefore,  may  be  attributed 
to  a variety  of  factors.  The  approach  taken  is  to  attribute  the  discrepancy 
to  the  most  obvious  physical  mechanism.  The  merit  of  this  assumption  can  then 
be  judged  by  the  final  results  obtained. 

Insight  may  be  gained  by  viewing  the  predicted  temperature  profiles.  Since 
radiative  and  conductive  heat  transfer  have  been  treated  separately,  the  temp- 
erature profiles  are  plotted  separately  with  the  caveat  that  since  about 
half  the  heat  transfer  is  due  to  radiation  and  half  is  due  to  conduction,  the 
initial  predicted  temperature  profile  may  be  roughly  taken  as  the  average  of 
the  two.  The  conduction  temperature  profile,  T]^ , is  taken  as  linear  from 
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T2  = 283.0  K to  Ti  = 310.3  K as  shown  in  Figure  6.  The  radiation  temperature 
profile  was  derived  by  Deissler  [40],  and  will  be  presented  here  in  brief. 

The  jump  in  emissivity  at  surface  2 is  given  as: 


®a2  = ^w2  - e„2)(l/£2  ~ 0.5)  (18) 

The  blackbody  emissivity  power  then  varies  linearly  (with  x)  through  the  medium, 
viz : 


e^(x)  = e^o  + X 0.75p3L  (^1-^2)  

^ 0.75p8L  + 1/e^  +1/S2-1 


(19) 


Combining  the  two  previous  equations,  the  temperature  profile  in  the  medium  is 
determined  as : 


Tr  = { 


"w2  ^®wl  - ^w2  ^ |~  1 _ Q 3 + 


0. 75pSL  X ^ 


0.75p3L  + 1/cj^  + l/e2 


(20) 


The  profiles  of  Tj-  are  shown  in  Figure  6 for  the  four  thicknesses  of  interest 
where  the  parameters  have  been  taken  as  = 9 kg»m  and  = C2  “ 0.9.  The 
jump  boundary  conditions  result  in  an  increase  above  T]^^  at  x = 0 and  a decrement 
below  T^  at  x =1.  There  is  clearly  a significant  increase  in  the  average  temper- 
ature when  radiation  is  considered,  and  it  is  observed  that  the  radiation  profiles 
are  nonlinear.  Conduction  will  have  the  following  two  effects  on  the  radiative 
temperature  profile:  (1)  it  will  eliminate  the  jump  at  the  boundaries,  and  (2) 

it  will  tend  to  smooth  the  overall  profile  by  transferring  heat  through  the  insu- 
lation. This  second  effect  will  serve  to  reduce  the  magnitude  of  the  temperature 
difference  from  conduction  and,  hence,  will  reduce  the  thickness  effect.  This 
is  precisely  the  type  of  mechanism  needed  to  enhance  the  present  theory. 
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Attention  is  therefore  focused  upon  the  coupling  between  the  radiative  and  the 
conductive  profiles.  Since  the  effects  are  small,  the  base  solutions  as  given 
above  will  be  used,  similar  to  perturbation  theory,  to  determine  the  coupled 
heat  transfer.  The  initial  approach  taken  was  to  integrate  the  radiative  temp- 
eratures excess  (as  a function  of  position)  divided  by  the  conduction  length  to 
the  cold  wall.  This  approach  did  not  bear  fruit  as  the  temperature  jump  at  the 
cold  wall  rendered  the  integrand  and  subsequently  the  integral  infinite.  Since 
the  conduction  will  eliminate  the  jump,  the  mathematics  is  not  inconsistent 
with  the  physics  of  the  situation.  Next,  it  was  reasoned  that  the  contribution 
to  the  effective  thermal  conductivity  of  the  coupling  between  radiation  and 
conduction  should  be  proportional  to  the  absolute  magnitude  of  the  average 
radiative  temperature  excess,  i.e.: 


x=l 

kcr  = T / 

x=0 


(25) 


This  integral,  denoted  as  the  non-dimensional  temperature  integral,  I'p,  is  shown 

in  Figure  7 , as  a function  of  density  for  the  four  thicknesses  of  interest  with 

2 —1 

= ^2  “ 3 = 4 m »kg  . The  graph  appears  to  be  consistent  with 

the  observed  data,  i.e.,  the  curves  decrease  with  increasing  density  and  the 
effects  are  more  pronounced  at  lower  thicknesses. 


The  approach  taken  is  to  apply  the  sequential  least  squares  filtering  technique 
to  estimate  y in  equation  (21).  Since  this  term  is  likely  to  be  small,  the 
filtering  is  done  for  this  term  alone  to  avoid  possible  computational  "condition- 
ing" problems.  Thus,  the  conduction  and  the  primary  radiation  terms  are  taken 
to  the  right-hand-side  of  the  equations  as  knowns , with  an  assumed  3.  The 
sequential  filtering  is  then  done  for  a variety  of  values  of  3,  resulting  in 
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a variance,  o^,  for  each  assumed  3.  This  corresponds  to  hunting  around 
in  the  two  space  of  (Yj3)  for  the  overall  minimum,  and  is  illustrated 
graphically  in  Figure  8.  The  minimum  is  achieved  for  8 = 44.0  m^'kg”!, 
and  the  corresponding  "best-fit"  value  of  Y if  Y = 0,0740  W*m”^»K“^. 

Recall  that  a correction  of  about  0.0015  W*m“l»K“^  is  being  sought, 
and  the  integral,  I'p,  was  shown  to  have  a value  of  about  0.025.  Thus, 

Y appears  to  be  of  the  correct  magnitude. 

The  new  relationship  is  presented  together  with  the  experimental  data  in 
Figure  9 for  the  four  thicknesses.  The  predicted  thickness  effect  is  now  seen 
to  be  consistent  with  the  experimental  data.  Moreover,  it  appears  that  the 
thickness  effect  is  delineated  about  as  much  as  possible  given  the  scatter  of 
the  data.  Two  datum  points  from  the  curve  fit  of  Rennex  [11]  at  the  two  thick- 
nesses of  0,0254  m (1")  and  0.1524  m (6")  are  shown  for  3 = 44  m2»kg~l  as 
open  and  closed  hexagons,  respectively.  The  thickness  effect  is  about  the 
same  as  is  predicted  by  the  present  analysis,  but  these  values  lie  well 
below  the  present  data.  The  lower  values  are  attributed  to  the  fact  that  the 
value  of  3 used  in  the  calculation  is  different  from  the  one  that  would  be 
obtained  by  filtering  using  his  curve  fit.  This  illustrates  that  the  "best-fit" 
parameters  for  one  model  are  not  necessarily  the  "best-fit"  values  for  another 
model.  Using  the  previous  value  of  3 = 39.5  m-*kg~^,  the  open  hexagons 
are  plotted  for  all  densities  and  at  the  thicknesses  of  0.0254  m (1")  and 
0.1524  m (6")  with  a double-dashed  line  and  a triple-dashed  line,  respectively. 
From  this,  it  is  observed  that  the  Rennex-curve  fit  slightly  overpredicts  the 
thickness  effect  for  the  present  data  set,  but  the  overall  agreement  is  good. 
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4.  CONCLUSIONS  AND  RECOMMENDATIONS 


Conduction  has  been  shown  to  be  essentially  that  of  stagnant  air,  and 
convection  has  been  demonstrated  to  be  nonexistent  in  horizontal  building 
insulation  under  the  usual  circumstances  encountered  in  practice.  Radiation 
was  found  to  contribute  approximately  one-half  to  the  total  heat  transfer.  A 
relation  of  the  diffusion  approximation  genre  was  found  to  be  inappropriate  in 
form  to  predict  the  thickness  effect.  The  inclusion  of  a term  to  account  for 
the  conductive-radiative  interaction  was  found  to  represent  the  data  well. 

The  present  approximation  was  found  to  be  slightly  better  for  the  present 
data  set  than  the  curve  fit  of  Rennex  [11].  For  completeness,  the  entire 
relationship  is  presented  as  follows: 

^c  = + ^r  + TIt  (22) 

where  the  conductive  contribution  is  given  as  [10]: 


^c  ~ ^g  *^k  ( ^ 


1 - 


1 + 


) 

2_Mv 


+ (1-cti^) 


(1 


1 - 


1 + 


Rr 

3(1  + R^) 


(23) 


and  the  radiative  term  is  as  follows  [40]: 


_ ^^^1^  ^ T2^)(T^  + T2)L 
^ ~ 3 P3L/4  + l/e^+l/£2-l 

The  values  of  the  constants  are  given  as  follows: 


3 = 

44.0  m^'kg  ^ 

(25) 

Y = 

0.0740  W*m“^*K~3 

(26) 

2^k  = 

0.5 

C7 
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The  present  approach  predicts  the  thickness  effect  very  well  from  a firm  fundamental 
basis.  The  universal  merit  of  the  present  result  can  only  be  assessed  when 
the  technique  has  been  applied  to  an  inclusive  data  set  including  different 
fiber  diameters  (most  especially),  different  surface  temperatures,  and  a wider 
variation  of  densities.  Also,  the  variation  in  types  and  quantities  of  binders 
is  a complicating  factor. 
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Table  1 . Experimental  Data 


Thickness 

(mm) 

(W.^Hr-I) 

Pm 

(kg*m“3) 

PT 

(kg*m  J) 

25.27 

0.04470 

10.03* 

10.03 

25.46 

0.04542 

9.47 

9.52 

25.35 

0.04618 

9.14 

8.92 

25.35 

0.04617 

9.11 

9.04 

25.35 

0.04557 

8.83 

8.78 

25.35 

0.04619 

9.25 

9.98 

25.35 

0.04810 

8.09 

8.67 

25.35 

0.04362 

10.67 

10.39 

25.35 

0.04443 

9.92 

9.91 

25.35 

0.04427 

9.17 

10.17 

25.38 

0.04487 

10.56* 

10.56 

25.35 

0.04512 

9.13 

9.09 

25.37 

0.04577 

9.11 

9.39 

25.35 

0.04568 

9.01 

9.22 

25.38 

0.04527 

9.63 

9.78 

25.38 

0.04457 

10.22 

10.79 

25.77 

0.04377 

11.13 

11.14 

25.77 

0.04364 

10.88 

11.07 

38.08 

0.04715 

9.51 

8.64 

38.08 

0.05030 

7.14 

7.34 

38.15 

0.04272 

11.40 

11.24 

76.26 

0.04949 

8.33 

8.42 

76.91 

0.04721 

9.20 

8.99 

76.91 

0.04968 

8.22 

8.34 

76.15 

0.04647 

9.27 

9.25 

76.15 

0.04777 

8.85 

8.56 

76.15 

0.04554 

9.75 

9.28 

76.15 

0.04736 

8.89 

9.06 

76.22 

0.04800 

8.78* 

8.78 

76.15 

0.04632 

9.51 

9.07 

76.22 

0.04717 

9.10 

9.18 

76.18 

0.04746 

8.95 

8.97 

152.35 

0.04752 

8.96 

9.06 

152.35 

0.04768 

8.94 

9.07 

152.35 

0.04798 

9.03 

9.00 

152.35 

0.04680 

9.36 

9.00 

152.35 

0.04880 

8.47 

8.92 

152.35 

0.04780 

8.95 

9.14 

152.35 

0.04788 

8.93 

9.16 

152.42 

0.04882 

8.94* 

8.94 

152.35 

0.04808 

8.78 

8.88 

152.37 

0.04772 

8.97 

8.88 

152.35 

0.04824 

8.69 

8.67 

152.42 

0.04713 

9.27 

9.32 

152.34 

0.04803 

8.93 

9.31 

152.35 

0.04752 

8.96 

9.06 

152.39 

0.04677 

9.16 

8.97 

* No  metering  sections  densities  were  available  for  these  points,  so 
the  metering  section  density  was  taken  as  the  total  specimen  density. 
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LEAST  SQUARES  DEVELOPMENT 


APPENDIX  A.  LEAST  SQUARES  DEVELOPMEIT 


A.l  INTRODUCTION 

A formal  mathematical  development  of  the  sequential  least  squares  algorithm  [43] 
used  in  this  work  is  presented  in  this  appendix.  In  the  following  discussion, 
matrices  will  be  devoted  by  upper  case  letters,  vectors  by  lower  case  letters 
with  a bar  above,  and  scalars  by  lower  case  letters  with  no  bar.  All  subscripts 
will  refer  to  order  in  time  or  in  sequence,  although  in  the  general  setting  of 
linear  systems  no  reference  to  time  is  necessary. 

A. 2 GENERAL  DEVELOPMENT 

The  general  over-determined  linear  system  is  given  by: 


where  A is  k x n (k>n)  and  x and  y are  of  dimension  n and  k respectively.  In 
general,  no  solution  exists  which  satisfies  equations  (A.l),  so  it  is  common 
practice  to  choose  as  a "solution"  that  value  of  x,  denoted  Xj^  g,  which 
minimizes  the  functional 


I.e.  j is  the  sum  of  the  squared  residuals. 

The  derivation  of  the  least  squares  solution  to  equation  (A.l)  is 
straightforward,  using  the  following  identities  for  an  arbitrary  matrix  B and 
vector  z: 


Ax  = y. 


(A.l) 


j = (Ax-y)T  (Ax-y) 


(A. 2) 


(A. 3) 


3(z^B)  ^ 

9z 


(A. 4) 


3(3z)  = bT 

9z 


A-1 


To  minimize  equation  (A. 2),  j is  differentiated  with  respect  to  x and  the 

9t 

derivative  ~ is  set  equal  to  zero, 

3x 


_ 9(x'^a'^Ax)  _ 9 (x’^a'^y)  _ 9(y^Ax)  _ 9 (y'^y) 

9x  ~5x  9x  9x  9x 

= 2x\'^Ax-  A^y  - (y^A)'^ 

= 2aTax-  2A^y 

= 0 ^ (A. 6) 

Therefore : 

x^^g  = (A^A)"^  A^y  (A. 7) 

A. 3 SEQUENTIAL  LEAST  SQUARES 

It  is  convenient  to  perform  the  operations  specified  by  equations  (A. 7)  in  a 
sequential  manner,  without  requiring  excessive  amounts  of  computer  storage  and 
without  necessitating  the  performance  of  a matrix  inverse  as  each  new  row  of 
the  system  given  by  equation  (A.l)  is  obtained.  To  this  end,  the  Woodbury 
matrix  inversion  lemma  is  introduced: 

(F+  BCD)"1  = F"1  - F^1b(C"1  + DF^1b)“1dF~1  (A, 8) 


It  is  now  useful  to  rewrite  equations  (A.l)  as  a block  system  as  follows: 


XI 

^2 

• 

X2 

• 

= 

^2 

• 

• 

• 

• 

Xn 

• 

^k 

; or  A^  X = y,^. 


(A. 9) 


A-2 


where  represents  the  i—  row  of  the  A matrix  and  represents  the  A matrix 
before  the  (k+1)^  row  has  been  adjoined  to  the  system. 


Assume  that  for  some  value  of  k,  the  system,  equation  (A. 9),  has  been  solved 
for  s ^^-en  the  information  and  become  available,  ^ can 

then  be  updated  to  ^ using  equation  (A. 8)  . 

First,  define  and  by: 

^k  " \ " ^k^k  (A.  10) 


Then  it  is  seen  that  for  any  value  of  k 
''^l.s,k  = ^k^k* 


(A. 11) 


To  update  and  to  P^+l  and  note  that 


1 _ 

aT 

3 T 

Ak 

^k+1 

\+l 

~ ^k  \ ^k+1  ^k+1 

" \ ^ ^k+1  ^k+1 

with  C = 1 in  equation  (A. 8),  = ^^k+1  found  as  follows: 

\+l  " \ ~ \ ^k+1*^^  ^k+1  ^k  ^k+1^  ^ ^k+1  ^k* 


(A. 12) 


(A. 13) 


Note  in  particular  that  the  quality  (1  + ai,^j^  P^  ^k+1^  ^ scaler,  so  that 

updating  Pj^  to  Pj^+i  requires  no  matrix  inversion. 


The  vector  z is  also  easily  updated  as  follows: 
^k+1  " ^k+1  ^k+1 


aT 

^k 

^k+1 

\+l 

A-3 


(A. 14) 


= Ak  7k  + 4+1  7k+l 

" 4+1  ^k+1 

Storage  requirements  for  the  above  sequential  least  squares  scheme  are  minimal. 
The  entire  A matrix  and  y vector  need  not  be  saved;  only  the  relatively  small  P 
matrix  (n  x n)  and  the  z vector  (dimensional  n)  need  to  be  saved.  As  new  data 
arrives  it  is  absorbed  by  the  current  P and  z,  and  the  least  squares  estimate 
of  X is  kept  current. 
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for  different  thicknesses  results  in  a parameter  estimate  of  the  optical 
extinction  coefficient.  This  parameter  determines  the  amount  of  change  in 
apparent  thermal  conductivity  as  a function  of  sample  thickness.  This 
phenomena  is  referred  to  as  the  "thickness  effect." 
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